library(dplyr)
library(ggplot2)
biopics <- read.csv("~/edi_imp/datacamp_JB/handing/biopics.csv")

Capítulo 1

Los datos faltantes están en todas partes. El proceso de rellenar los valores faltantes se conoce como imputación, y saber cómo rellenar correctamente los datos faltantes es una habilidad esencial si desea producir predicciones precisas y distinguirse de la multitud. En este curso, aprenderá a usar visualizaciones y pruebas estadísticas para reconocer patrones de datos faltantes y cómo imputar datos utilizando una colección de modelos estadísticos y de aprendizaje automático. También adquirirá habilidades para la toma de decisiones, lo que le ayudará a decidir qué método de imputación se adapta mejor a una situación particular. Por último, aprenderá a incorporar la incertidumbre de la imputación en su inferencia y predicciones, haciéndolos más robustos y confiables.

Sesión 1

Regresión lineal con datos incompletos

Los datos faltantes son un problema común y tratarlos de manera adecuada es extremadamente importante. Ignorar los puntos de datos faltantes o completarlos incorrectamente puede hacer que los modelos funcionen de maneras inesperadas y causar que las predicciones e inferencias estén sesgadas.

En este capítulo, trabajará con el conjunto de datos biopics. Contiene información sobre varias películas biográficas, incluyendo sus ingresos, características del tema y algunas otras variables. Sin embargo, algunos de los puntos de datos están faltando. Los datos originales vienen con el paquete R fivethirtyeight, pero en este curso, trabajará con una versión ligeramente preprocesada.

En este ejercicio, conocerá el conjunto de datos y ajustará un modelo de regresión lineal para explicar los ingresos de una película. ¡Comencemos!

Actividad

Muestra las primeras 10 observaciones de los datos biopics y familiarízate con las variables.

# Print first 10 observations
head(biopics, 10)
##    country year earnings sub_num sub_type sub_race non_white sub_sex
## 1       UK 1971       NA       1 Criminal     <NA>         0    Male
## 2    US/UK 2013   56.700       1    Other  African         1    Male
## 3    US/UK 2010   18.300       1  Athlete     <NA>         0    Male
## 4   Canada 2014       NA       1    Other    White         0    Male
## 5       US 1998    0.537       1    Other     <NA>         0    Male
## 6       US 2008   81.200       1    Other    other         1    Male
## 7       UK 2002    1.130       1 Musician    White         0    Male
## 8       US 2013   95.000       1  Athlete  African         1    Male
## 9       US 1994   19.600       1  Athlete     <NA>         0    Male
## 10   US/UK 1987    1.080       2   Author     <NA>         0    Male
# Get the number of missing values per variable
biopics %>%
    is.na() %>% 
    colSums()
##   country      year  earnings   sub_num  sub_type  sub_race non_white   sub_sex 
##         0         0       324         0         0       197         0         0
# Fit linear regression to predict earnings
model_1 <- lm(earnings ~ country + year + sub_type, 
              data = biopics)

# Fit linear regression to predict earnings
model_2 <- lm(earnings ~ country + year + sub_type + sub_race, 
              data = biopics)

Reconociendo los mecanismos de datos faltantes

En este ejercicio, se presentarán seis escenarios diferentes en los que faltan algunos datos. Intenta asignar a cada uno de ellos el mecanismo de datos faltantes más probable. Como recordatorio, aquí hay algunas pautas generales:

Si la razón de la falta de datos es puramente aleatoria, es MCAR. Si la razón de la falta de datos puede explicarse por otra variable, es MAR. Si la razón de la falta de datos depende del valor faltante en sí mismo, es MNAR.

alt text

alt text

Prueba t para MAR: preparación de datos

¡Buen trabajo al clasificar los mecanismos de datos faltantes en el último ejercicio! De los tres, MAR es posiblemente el más importante de detectar, ya que muchos métodos de imputación asumen que los datos son MAR. Este ejercicio, por lo tanto, se centrará en la prueba de MAR.

Trabajarás con los datos familiares de biopics. El objetivo es probar si el número de valores faltantes en earnings difiere por género del sujeto. En este ejercicio, solo prepararás los datos para la prueba t. Primero, crearás una variable ficticia que indique la falta de datos en earnings. Luego, la dividirás por género filtrando los datos para mantener uno de los géneros y luego sacando la variable ficticia. Para filtrar, puede ser útil imprimir el head() de biopics en la consola y examinar la variable de género.

# Create a dummy variable for missing earnings
biopics <- biopics %>% 
  mutate(missing_earnings = is.na(earnings))

# Pull the missing earnings dummy for males
missing_earnings_males <- biopics %>% 
  filter(sub_sex == "Male") %>% 
  pull(missing_earnings)

# Pull the missing earnings dummy for females
missing_earnings_females <- biopics %>% 
  filter(sub_sex == "Female") %>% 
  pull(missing_earnings)

Prueba t para MAR: interpretación

En el ejercicio anterior, has preparado dos vectores con los valores faltantes de ingresos para cada género: missing_earnings_males y missing_earnings_females. Ambos están disponibles en tu espacio de trabajo. ¡Ahora puedes realizar la prueba t para verificar si sus medias difieren significativamente entre sí! ¡Realicemos algunas pruebas estadísticas serias!

# Run the t-test
t.test(missing_earnings_males, missing_earnings_females)
## 
##  Welch Two Sample t-test
## 
## data:  missing_earnings_males and missing_earnings_females
## t = 1.1116, df = 294.39, p-value = 0.2672
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.03606549  0.12969214
## sample estimates:
## mean of x mean of y 
## 0.4366438 0.3898305

Sesión 2

Aggregation plot

El gráfico de agregación proporciona la respuesta a la pregunta básica que uno puede hacer sobre un conjunto de datos incompleto: ¿en qué combinaciones de variables faltan datos y con qué frecuencia? Es muy útil para obtener una visión general de alto nivel de los patrones de ausencia de datos. Por ejemplo, hace visible inmediatamente si hay alguna combinación de variables que faltan juntas con frecuencia, lo que podría sugerir alguna relación entre ellas.

En este ejercicio, primero dibujarás el gráfico de agregación para los datos de biopics y luego practicarás sacando conclusiones basadas en él. ¡Vamos a hacer algunos gráficos!

# Load the VIM package
library(VIM)
## Loading required package: colorspace
## Loading required package: grid
## VIM is ready to use.
## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues
## 
## Attaching package: 'VIM'
## The following object is masked from 'package:datasets':
## 
##     sleep
# Draw an aggregation plot of biopics
biopics %>% 
    aggr(combined = TRUE, numbers = TRUE)

Preguntas

Basado en el gráfico de agregación que acaba de crear, ¿cuál de las siguientes afirmaciones es falsa?

Posibles respuestas:

El 10% de las observaciones tienen valores faltantes tanto en earnings como en sub_race.

Hay más valores faltantes en sub_race que en earnings.

El 42% de las observaciones no tiene entradas faltantes.

Exactamente dos variables en los datos biopics tienen valores faltantes.

Sesión 3

gráfico de columna vertebral

El gráfico de agregación que has dibujado en el ejercicio anterior te dio una visión general de alto nivel de los datos faltantes. Si estás interesado en la interacción entre variables específicas, un gráfico de columna vertebral es el camino a seguir. Te permite estudiar el porcentaje de valores faltantes en una variable para diferentes valores de la otra, lo cual es conceptualmente muy similar a los test t que has estado realizando en la lección anterior.

En este ejercicio, dibujarás un gráfico de columna vertebral para investigar el porcentaje de datos faltantes en earnings para diferentes categorías de sub_race. ¿Hay más datos faltantes en earnings para algunas razas específicas del personaje principal de la película? ¡Vamos a descubrirlo! El paquete VIM ya ha sido cargado para ti.

# Draw a spine plot to analyse missing values in earnings by sub_race

biopics %>% 
    dplyr::select(sub_race,earnings) %>%
    spineMiss()

### Preguntas

Basándose en la gráfica de columna vertebral que acabas de crear, ¿cuál de las siguientes afirmaciones es falsa?

Opciones de respuesta:

En la gran mayoría de las películas, el personaje principal es blanco.

Cuando el sujeto principal es africano, es más probable que tengamos información completa sobre las ganancias.

En lo que respecta a las ganancias y la subraza, los datos parecen ser MAR.

La raza que aparece con menos frecuencia en los datos tiene alrededor del 40% de las ganancias faltantes. (incorrecta)

Mosaic plot

La gráfica de columna vertebral que has creado en el ejercicio anterior te permite estudiar los patrones de datos faltantes entre dos variables a la vez. Esta idea se generaliza a más variables en forma de un gráfico de mosaico.

En este ejercicio, comenzarás por crear una variable ficticia que indique si Estados Unidos participó en la producción de cada película. Para hacer esto, utilizarás la función grepl(), que verifica si la cadena pasada como su primer argumento está presente en el objeto pasado como su segundo argumento. Luego, dibujarás un gráfico de mosaico para ver si el género del sujeto se correlaciona con la cantidad de datos faltantes en earnings tanto para películas estadounidenses como no estadounidenses.

Los datos biopics y el paquete VIM ya están cargados para ti. ¡Vamos a hacer algunas gráficas exploratorias!

Ten en cuenta que se ha creado una función de propiedad display_image() para devolver la salida de la última versión del paquete VIM. Asegúrate de expandir la sección Visor HTML.

# Prepare data for plotting and draw a mosaic plot
biopics %>%
    # Create a dummy variable for US-produced movies
    mutate(is_US_movie = grepl("US", country)) %>%
    # Draw mosaic plot
    mosaicMiss(highlight = "earnings", 
             plotvars = c("is_US_movie", "sub_sex"))

# Return plot from latest VIM package - expand the HTML viewer section
#display_image()

Capítulo 2

Sesión 1

Olfateando el peligro de la imputación media

Uno de los métodos de imputación más populares es la imputación media, en la cual los valores faltantes en una variable se reemplazan con la media de los valores observados en esta variable. Sin embargo, en muchos casos, este enfoque simple es una mala elección. A veces, una mirada rápida a los datos puede alertarte sobre los peligros de la imputación media.

En este capítulo, trabajarás con una submuestra de los datos del proyecto de Atmósfera Tropical Oceánica (tao). El conjunto de datos consiste en mediciones atmosféricas tomadas en dos períodos de tiempo diferentes en cinco ubicaciones distintas. Los datos vienen con el paquete VIM.

En este ejercicio, te familiarizarás con los datos y realizarás un análisis simple que indicará cuáles podrían ser las consecuencias de la imputación media. ¡Echemos un vistazo a los datos de tao!

data(tao, package = "VIM")
names(tao)<-tolower(names(tao))
names(tao)<-sub("[.]", "_", names(tao))
names(tao)<-sub("[.]", "_", names(tao))

# Print first 10 observations
head(tao, 10)
##    year latitude longitude sea_surface_temp air_temp humidity uwind vwind
## 1  1997        0      -110            27.59    27.15     79.6  -6.4   5.4
## 2  1997        0      -110            27.55    27.02     75.8  -5.3   5.3
## 3  1997        0      -110            27.57    27.00     76.5  -5.1   4.5
## 4  1997        0      -110            27.62    26.93     76.2  -4.9   2.5
## 5  1997        0      -110            27.65    26.84     76.4  -3.5   4.1
## 6  1997        0      -110            27.83    26.94     76.7  -4.4   1.6
## 7  1997        0      -110            28.01    27.04     76.5  -2.0   3.5
## 8  1997        0      -110            28.04    27.11     78.3  -3.7   4.5
## 9  1997        0      -110            28.02    27.21     78.6  -4.2   5.0
## 10 1997        0      -110            28.05    27.25     76.9  -3.6   3.5
# Get the number of missing values per column
tao %>%
  is.na() %>% 
  colSums()
##             year         latitude        longitude sea_surface_temp 
##                0                0                0                3 
##         air_temp         humidity            uwind            vwind 
##               81               93                0                0
# Calculate the number of missing values in air_temp per year
tao %>% 
  group_by(year) %>% 
  summarize(num_miss = sum(is.na(air_temp)))
## # A tibble: 2 × 2
##    year num_miss
##   <int>    <int>
## 1  1993        4
## 2  1997       77

Imputación de la media en temperatura

Imputar la media en la temperatura puede ser arriesgado. Si la variable que se está imputando está correlacionada con otras variables, esta correlación podría ser destruida por los valores imputados. Lo viste en el ejercicio anterior cuando analizaste la variable air_temp.

Para averiguar si estas preocupaciones son válidas, en este ejercicio realizarás una imputación de la media en air_temp, creando también un indicador binario para mostrar dónde se imputan los valores. Será útil en el siguiente ejercicio, cuando evaluarás el desempeño de tu imputación. ¡Vamos a completar esos valores faltantes!

tao_imp <- tao %>% 
  # Create a binary indicator for missing values in air_temp
  mutate(air_temp_imp = ifelse(is.na(air_temp), TRUE, FALSE)) %>%
  # Impute air_temp with its mean
  mutate(air_temp = ifelse(is.na(air_temp), mean(air_temp, na.rm = TRUE), air_temp))

# Print the first 10 rows of tao_imp
head(tao_imp, 10)
##    year latitude longitude sea_surface_temp air_temp humidity uwind vwind
## 1  1997        0      -110            27.59    27.15     79.6  -6.4   5.4
## 2  1997        0      -110            27.55    27.02     75.8  -5.3   5.3
## 3  1997        0      -110            27.57    27.00     76.5  -5.1   4.5
## 4  1997        0      -110            27.62    26.93     76.2  -4.9   2.5
## 5  1997        0      -110            27.65    26.84     76.4  -3.5   4.1
## 6  1997        0      -110            27.83    26.94     76.7  -4.4   1.6
## 7  1997        0      -110            28.01    27.04     76.5  -2.0   3.5
## 8  1997        0      -110            28.04    27.11     78.3  -3.7   4.5
## 9  1997        0      -110            28.02    27.21     78.6  -4.2   5.0
## 10 1997        0      -110            28.05    27.25     76.9  -3.6   3.5
##    air_temp_imp
## 1         FALSE
## 2         FALSE
## 3         FALSE
## 4         FALSE
## 5         FALSE
## 6         FALSE
## 7         FALSE
## 8         FALSE
## 9         FALSE
## 10        FALSE

Evaluar la calidad de la imputación con un gráfico de margen

En el último ejercicio, has imputado la media de air_temp y has agregado una variable indicadora para denotar cuáles valores fueron imputados, llamada air_temp_imp. Ahora es momento de ver qué tan bien funciona esto.

Al examinar los datos de tao, podrías haber notado que también contiene una variable llamada sea_surface_temp, que razonablemente se esperaría que esté positivamente correlacionada con air_temp. Si ese es el caso, esperarías que estas dos temperaturas sean altas o bajas al mismo tiempo. Imputar la temperatura media del aire cuando la temperatura del mar es alta o baja rompería esta relación.

Para averiguarlo, en este ejercicio seleccionarás las dos variables de temperatura y la variable indicadora y las usarás para dibujar un gráfico de margen. ¡Vamos a evaluar la imputación media!

# Draw a margin plot of air_temp vs sea_surface_temp
tao_imp %>% 
  select(air_temp, sea_surface_temp, air_temp_imp) %>%
  marginplot(delimiter = "imp")

Sesión 2

Vanilla hot-deck

La imputación por hot-deck es un método simple que reemplaza cada valor faltante en una variable por el último valor observado en esa variable. Es muy rápido, ya que solo se necesita una pasada por los datos, pero en su forma más simple, hot-deck a veces puede romper las relaciones entre las variables.

En este ejercicio, lo probarás en el conjunto de datos tao. Imputarás los valores faltantes en la columna de temperatura del aire air_temp por hot-deck y luego dibujarás un gráfico de margen para analizar la relación entre los valores imputados y la columna de temperatura de la superficie del mar sea_surface_temp. ¡Veamos cómo funciona!

# Load VIM package
library(VIM)

# Impute air_temp in tao with hot-deck imputation
tao_imp <- hotdeck(tao, variable = "air_temp")

# Check the number of missing values in each variable
tao_imp %>% 
    is.na() %>% 
    colSums()
##             year         latitude        longitude sea_surface_temp 
##                0                0                0                3 
##         air_temp         humidity            uwind            vwind 
##                0               93                0                0 
##     air_temp_imp 
##                0
# Draw a margin plot of air_temp vs sea_surface_temp
tao_imp %>% 
    select(air_temp, sea_surface_temp, air_temp_imp) %>% 
    marginplot(delimiter = "imp")

¿Se ve bien la imputación? Observa las observaciones en la parte superior izquierda del gráfico con air_temp imputados y alta sea_surface_temp. Estas observaciones deben haber sido precedidas por observaciones con baja air_temp en el marco de datos, y por lo tanto, después de la imputación hot-deck, terminaron siendo valores atípicos con baja air_temp y alta sea_surface_temp.

Hot-deck trucos y consejos I: imputando dentro de dominios

Un truco que puede ayudar cuando la imputación por hot-deck rompe las relaciones entre las variables es imputar dentro de dominios. Esto significa que si la variable a imputar está correlacionada con otra variable categórica, se puede ejecutar hot-deck por separado para cada una de sus categorías.

Por ejemplo, se podría esperar que la temperatura del aire dependa del tiempo, ya que estamos viendo que las temperaturas promedio aumentan debido al calentamiento global. El indicador de tiempo que tienes disponible en los datos de tao es una variable categórica, año. Primero, comprobaremos si la temperatura media del aire es diferente en cada uno de los dos años estudiados y luego ejecutaremos hot-deck dentro de los dominios de los años. Finalmente, volveremos a dibujar el gráfico de márgenes para evaluar el rendimiento de la imputación.

# Calculate mean air_temp per year
tao %>% 
    group_by(year) %>% 
    summarize(average_air_temp = mean(air_temp, na.rm = TRUE))
## # A tibble: 2 × 2
##    year average_air_temp
##   <int>            <dbl>
## 1  1993             23.4
## 2  1997             27.1
# Hot-deck-impute air_temp in tao by year domain
tao_imp <- hotdeck(tao, variable = "air_temp", domain_var = "year")

# Draw a margin plot of air_temp vs sea_surface_temp
tao_imp %>% 
    select(air_temp, sea_surface_temp, air_temp_imp) %>% 
    marginplot(delimiter = "imp")

Los resultados se ven mucho mejor esta vez. Sin embargo, si observas la esquina superior derecha del gráfico, verás que la varianza en los valores imputados (naranja) es algo mayor que entre los valores observados (azul). ¡Veamos si podemos mejorar aún más en el próximo ejercicio!

Hot-deck trucos y consejos II: ordenando por variables correlacionadas

Otro truco que puede mejorar el rendimiento de la imputación hot-deck es ordenar los datos por variables correlacionadas con la que queremos imputar.

Por ejemplo, en todos los gráficos de márgenes que ha estado dibujando recientemente, ha visto que la temperatura del aire está fuertemente correlacionada con la temperatura de la superficie del mar, lo cual tiene mucho sentido. Puede aprovechar este conocimiento para mejorar su imputación hot-deck. Si primero ordena los datos por sea_surface_temp, entonces cada valor imputado de air_temp vendrá de un donante con una sea_surface_temp similar. ¡Veamos cómo funcionará esto!

# Hot-deck-impute air_temp in tao ordering by sea_surface_temp
tao_imp <- hotdeck(tao, variable = "air_temp", ord_var = "sea_surface_temp")

# Draw a margin plot of air_temp vs sea_surface_temp
tao_imp %>% 
    select(air_temp, sea_surface_temp, air_temp_imp) %>% 
    marginplot(delimiter = "imp")

Esta vez la imputación parece no afectar la relación entre las temperaturas del aire y la superficie del mar: si no fuera por los colores, probablemente no sabrías cuáles son los valores imputados. La imputación hot-deck, posiblemente mejorada con la imputación por dominios o el ordenamiento, es un método rápido y sencillo que puede funcionar bien en muchas situaciones. Sin embargo, a veces puede ser necesario un enfoque más complejo.

Sesión 3

Elegir el número de vecinos

La imputación de k-Nearest-Neighbors (o kNN) completa los valores faltantes en una observación en función de los valores que provienen de las k otras observaciones más similares a ella. El número de estas observaciones similares, llamadas vecinos, que se consideran es un parámetro que debe elegirse de antemano.

¿Cómo elegir k? Una forma es probar diferentes valores y ver cómo afectan las relaciones entre los datos imputados y observados.

Intentemos completar humidity en los datos de tao utilizando tres valores diferentes de k y ver cómo se ajustan los valores imputados a la relación entre humidity y sea_surface_temp.

Completar humedad con la imputación de kNN usando 30 vecinos y dibujar un marginplot() de sea_surface_temp vs humidity.

# Impute humidity using 30 neighbors
tao_imp <- kNN(tao, k = 30, variable = "humidity")

# Draw a margin plot of sea_surface_temp vs humidity
tao_imp %>% 
    select(sea_surface_temp, humidity, humidity_imp) %>% 
    marginplot(delimiter = "imp", main = "k = 30")

Imputa humidity con imputación kNN usando 15 vecinos y dibujando marginplot de sea_surface_temp vs humidity.

# Impute humidity using 15 neighbors
tao_imp <- kNN(tao, k = 15, variable = "humidity")

# Draw a margin plot of sea_surface_temp vs humidity
tao_imp %>% 
    select(sea_surface_temp, humidity, humidity_imp) %>% 
    marginplot(delimiter = "imp", main = "k = 15")

Imputa humidity con imputación kNN usando 5 vecinos y dibujando marginplot de sea_surface_temp vs humidity.

# Impute humidity using 5 neighbors
tao_imp <- kNN(tao, k = 5, variable = "humidity")

# Draw a margin plot of sea_surface_temp vs humidity
tao_imp %>% 
    select(sea_surface_temp, humidity, humidity_imp) %>% 
    marginplot(delimiter = "imp", main = "k = 5")

kNN trucos y consejos I: ponderando los donantes

Una variación de la imputación kNN que se aplica con frecuencia utiliza la llamada agregación ponderada por distancia. Lo que esto significa es que cuando agregamos los valores de los vecinos para obtener un reemplazo para un valor faltante, lo hacemos usando la media ponderada y las ponderaciones son las distancias invertidas de cada vecino. Como resultado, los vecinos más cercanos tienen más impacto en el valor imputado.

En este ejercicio, aplicarás la agregación ponderada por distancia mientras imputas los datos de tao. Esto solo requerirá pasar dos argumentos adicionales a la función kNN(). ¡Probémoslo!

# Load the VIM package
library(VIM)

# Impute humidity with kNN using distance-weighted mean
tao_imp <- kNN(tao, 
               k = 5, 
               variable = "humidity", 
               numFun = weighted.mean,
               weightDist = TRUE)

tao_imp %>% 
    select(sea_surface_temp, humidity, humidity_imp) %>% 
    marginplot(delimiter = "imp")

### Trucos y consejos de kNN II: ordenar variables

Mientras el algoritmo de k-Nearest Neighbors recorre las variables en los datos para imputarlos, calcula las distancias entre observaciones utilizando otras variables, algunas de las cuales ya han sido imputadas en los pasos anteriores. Esto significa que si las variables ubicadas al principio de los datos tienen muchas valores faltantes, entonces el cálculo de la distancia posterior se basa en muchos valores imputados. Esto introduce ruido en el cálculo de la distancia.

Por esta razón, es una buena práctica ordenar las variables cada vez más por el número de valores faltantes antes de realizar la imputación kNN. De esta manera, cada cálculo de distancia se basa en tantos datos observados y tan pocos datos imputados como sea posible.

¡Vamos a probar esto en los datos de tao!

# Get tao variable names sorted by number of NAs
vars_by_NAs <- tao %>%
  is.na() %>%
  colSums() %>%
  sort(decreasing = FALSE) %>% 
  names()

# Sort tao variables and feed it to kNN imputation
tao_imp <- tao %>% 
  select(vars_by_NAs) %>% 
  kNN(k= 5)
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(vars_by_NAs)` instead of `vars_by_NAs` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
tao_imp %>% 
    select(sea_surface_temp, humidity, humidity_imp) %>% 
    marginplot(delimiter = "imp")

El kNN que acabas de programar debería ser más preciso y resistente a imputaciones defectuosas, así que recuerda ordenar tus variables primero antes de realizar la imputación con kNN. ¡Esto nos lleva al final de este capítulo! ¡Sigue adelante! ¡Nos vemos en el Capítulo 3, donde aprenderás a usar modelos estadísticos y de aprendizaje automático para imputar valores faltantes!

Capítulo 3

Sesión 1

Imputación con regresión lineal

A veces, se puede utilizar el conocimiento del dominio, la investigación previa o simplemente el sentido común para describir las relaciones entre las variables en sus datos. En tales casos, la imputación basada en modelos es una gran solución, ya que le permite imputar cada variable de acuerdo con un modelo estadístico que puede especificar usted mismo, teniendo en cuenta cualquier suposición que pueda tener sobre cómo las variables impactan entre sí.

Para variables continuas, una elección de modelo popular es la regresión lineal. ¡No te restringe a relaciones lineales! Siempre puede incluir un cuadrado o un logaritmo de una variable en los predictores. En este ejercicio, trabajará con el paquete simputation para ejecutar una sola imputación de regresión lineal en los datos tao y analizar los resultados. ¡Vamos a intentarlo!

# Load the simputation package
library(simputation)

# Impute air_temp and humidity with linear regression
formula <- air_temp + humidity ~ year + latitude + sea_surface_temp
tao_imp <- impute_lm(tao, formula)
# Load the simputation package
library(simputation)

# Impute air_temp and humidity with linear regression
formula <- air_temp + humidity ~ year + latitude + sea_surface_temp
tao_imp <- impute_lm(tao, formula)

# Check the number of missing values per column
tao_imp %>% 
  is.na() %>% 
  colSums()
##             year         latitude        longitude sea_surface_temp 
##                0                0                0                3 
##         air_temp         humidity            uwind            vwind 
##                3                2                0                0
# Print rows of tao_imp in which air_temp or humidity are still missing 
tao_imp %>% 
  filter(is.na(air_temp) | is.na(humidity))
##   year latitude longitude sea_surface_temp air_temp humidity uwind vwind
## 1 1993        0       -95               NA       NA       NA  -5.6   3.1
## 2 1993        0       -95               NA       NA       NA  -6.3   0.5
## 3 1993       -2       -95               NA       NA     89.9  -3.4   2.4

La regresión lineal falla cuando al menos uno de los predictores está ausente. En este caso, fue sea_surface_temp. En el próximo ejercicio, lo solucionarás inicializando los valores faltantes antes de ejecutar impute_lm().

Inicialización de valores perdidos e iteración sobre variables

Como acabas de ver, la ejecución de impute_lm() podría no llenar todos los valores perdidos. Para asegurarte de imputar todos ellos, deberías inicializar los valores perdidos con un método simple, como la imputación de hot-deck que aprendiste en el capítulo anterior, que simplemente retroalimenta el último valor observado.

Además, una sola imputación generalmente no es suficiente. Se basa en los valores iniciales básicos y podría estar sesgada. Un enfoque adecuado es iterar sobre las variables, imputándolas una a la vez en las ubicaciones donde originalmente faltaban.

En este ejercicio, primero inicializarás los valores perdidos con la imputación de hot-deck y luego iterarás cinco veces sobre air_temp y humidity de los datos tao para imputarlos con regresión lineal. ¡Manos a la obra!

# Initialize missing values with hot-deck
tao_imp <- hotdeck(tao)

# Create boolean masks for where air_temp and humidity are missing
missing_air_temp <- tao_imp$air_temp_imp
missing_humidity <- tao_imp$humidity_imp

for (i in 1:5) {
  # Set air_temp to NA in places where it was originally missing and re-impute it
  tao_imp$air_temp[missing_air_temp] <- NA
  tao_imp <- impute_lm(tao_imp, air_temp ~ year + latitude + sea_surface_temp + humidity)
  # Set humidity to NA in places where it was originally missing and re-impute it
  tao_imp$humidity[missing_humidity] <- NA
  tao_imp <- impute_lm(tao_imp, humidity ~ year + latitude + sea_surface_temp + air_temp)
}

¡Esa es una aproximación profesional a la imputación basada en modelos que acabas de codificar! Pero, ¿cómo sabemos que 5 es el número adecuado de iteraciones para ejecutar? ¡Veamos la convergencia en el siguiente ejercicio!

Detectando convergencia

¡Excelente trabajo iterando sobre las variables en el ejercicio anterior! ¿Pero cuántas iteraciones son necesarias? Cuando los valores imputados no cambian con la nueva iteración, podemos detenernos.

Ahora extenderás tu código para calcular las diferencias entre las variables imputadas en las iteraciones subsiguientes. Para hacer esto, usarás la función de cambio porcentual absoluto promedio, definida para ti de la siguiente manera:

mapc <- function(a, b) { mean(abs(b - a) / a, na.rm = TRUE) }

mapc() es una función que te devuelve un solo número que te dice cuánto difiere b de a. La usarás para verificar cuánto cambian las variables imputadas en las iteraciones siguientes. En base a esto, decidirás cuántas iteraciones son necesarias.

Las máscaras booleanas missing_air_temp y missing_humidity están disponibles para ti, al igual que los datos de tao_imp inicializados con hot-deck.

mapc<- function(a, b) {
  mean(abs(b - a) / a, na.rm = TRUE)
}
diff_air_temp <- c()
diff_humidity <- c()

for (i in 1:5) {
  # Assign the outcome of the previous iteration (or initialization) to prev_iter
  prev_iter <- tao_imp
  # Impute air_temp and humidity at originally missing locations
  tao_imp$air_temp[missing_air_temp] <- NA
  tao_imp <- impute_lm(tao_imp, air_temp ~ year + latitude + sea_surface_temp + humidity)
  tao_imp$humidity[missing_humidity] <- NA
  tao_imp <- impute_lm(tao_imp, humidity ~ year + latitude + sea_surface_temp + air_temp)
  # Calculate MAPC for air_temp and humidity and append them to previous iteration's MAPCs
  diff_air_temp <- c(diff_air_temp, mapc(prev_iter$air_temp, tao_imp$air_temp))
  diff_humidity <- c(diff_humidity, mapc(prev_iter$humidity, tao_imp$humidity))
}

Question

Based on the differences stored in diff_air_temp and diff_humidity, what is the sufficient number of iterations to run?

To answer this question, you can print the two vectors in the console and analyze the numbers, or plot them using the function provided for you: just run plot_diffs(diff_air_temp, diff_humidity) in the console.

plot_diffs <- function(a, b) {
  data.frame("mapc" = c(a, b),
             "Variable" = c(rep("air_temp", length(a)),
                            rep("humidity", length(b))),
             "Iterations" = c(1:length(a), 1:length(b))) %>% 
    ggplot(aes(Iterations, mapc, color = Variable)) +
    geom_line(size = 1.5) +
    ylab("Mean absolute percentage change") +
    ggtitle("Changes in imputed variables' values across iterations") +
    theme(legend.position = "bottom")
}
plot_diffs(diff_air_temp, diff_humidity)

## Sesión 2

###imputación de regresión logística

Una opción popular para imputar variables binarias es la regresión logística. Desafortunadamente, no hay una función similar a impute_lm() que lo haga. ¡Por eso escribirás una función tú mismo!

Llamemos a la función impute_logreg(). Su primer argumento será un marco de datos df, cuyos valores faltantes se han inicializado y solo contiene valores faltantes en la columna a imputar. El segundo argumento será una fórmula para el modelo de regresión logística.

La función hará lo siguiente:

Mantendrá las ubicaciones de los valores faltantes.

Construirá el modelo.

Realizará predicciones.

Reemplazará los valores faltantes con las predicciones.

No te preocupes por la línea que crea imp_var - esto es solo una forma de extraer el nombre de la columna a imputar de la fórmula. ¡Vamos a hacer algo de programación funcional!

impute_logreg <- function(df, formula) {
  # Extract name of response variable
  imp_var <- as.character(formula[2])
  # Save locations where the response is missing
  missing_imp_var <- is.na(df[imp_var])
  # Fit logistic regression mode
  logreg_model <- glm(formula, data = df, family = binomial)
  # Predict the response and convert it to 0s and 1s
  preds <- predict(logreg_model, type = "response")
  preds <- ifelse(preds >= 0.5, 1, 0)
  # Impute missing values with predictions
  df[missing_imp_var, imp_var] <- preds[missing_imp_var]
  return(df)
}

La función que escribiste está completamente operativa y se puede enchufar en el bucle sobre las variables que viste en el capítulo anterior, al igual que impute_lm() del paquete simputation. Pronto, combinarás estos dos para imputar tanto variables continuas como binarias. Pero antes, mejoraremos tu impute_logreg() para que reproduzca mejor la variabilidad en los datos imputados.

Dibujando desde una distribución condicional

Simplemente llamar a predict() en un modelo siempre devolverá el mismo valor para los mismos valores de los predictores. Esto da como resultado una pequeña variabilidad en los datos imputados. Para aumentarla y que la imputación replique la variabilidad de los datos originales, podemos extraer de la distribución condicional. Esto significa que en lugar de siempre predecir 1 cuando el modelo devuelve una probabilidad mayor que 0,5, podemos extraer la predicción de una distribución binomial descrita por la probabilidad devuelta por el modelo.

Trabajarás en el código que escribiste en el ejercicio anterior. La siguiente línea fue eliminada:

preds <- ifelse(preds >= 0.5, 1, 0)

Tu tarea es llenar su lugar con la creación de una distribución binomial. ¡Eso es solo una línea de código!

 impute_logreg <- function(df, formula) {
  # Extract name of response variable
  imp_var <- as.character(formula[2])
  # Save locations where the response is missing
  missing_imp_var <- is.na(df[imp_var])
  # Fit logistic regression mode
  logreg_model <- glm(formula, data = df, family = binomial)
  # Predict the response
  preds <- predict(logreg_model, type = "response")
  # Sample the predictions from binomial distribution
  # preds <- ifelse(preds >= 0.5, 1, 0)
  preds <- rbinom(length(preds), size = 1, prob = preds)
  # Impute missing values with predictions
  df[missing_imp_var, imp_var] <- preds[missing_imp_var]
  return(df)
}

Dibujar la distribución condicional hará que la variabilidad de los datos imputados sea más similar a la del conjunto de datos observados originales. Con esta potente función en tus manos, ahora estás listo para diseñar un flujo de imputación basado en modelos que se encargue tanto de variables continuas como binarias. ¡Vamos a hacerlo en el siguiente ejercicio!

Imputación basada en modelos con varios tipos de variables

¡Buen trabajo al escribir la función para implementar la imputación de regresión logística con dibujo de la distribución condicional! ¡Has codificado estadísticas bastante avanzadas! En este ejercicio, combinarás lo que has aprendido hasta ahora sobre imputación basada en modelos para imputar diferentes tipos de variables en los datos de tao.

Tu tarea es iterar sobre las variables como lo has hecho en el capítulo anterior e imputar dos variables:

is_hot, una nueva variable binaria que se creó a partir de air_temp, que es 1 si air_temp está a 26 grados o más y 0 de lo contrario; humidity, una variable continua con la que ya estás familiarizado. Tendrás que utilizar la función de regresión lineal que aprendiste antes, así como tu propia función para la regresión logística. ¡Vamos a ello!

tao$is_hot<-ifelse(tao$air_temp>= 26, 1,0)
# Initialize missing values with hot-deck
tao_imp <- hotdeck(tao)

# Create boolean masks for where is_hot and humidity are missing
missing_is_hot <- tao_imp$is_hot_imp
missing_humidity <- tao_imp$humidity_imp

for (i in 1:3) {
  # Set is_hot to NA in places where it was originally missing and re-impute it
  tao_imp$is_hot[missing_is_hot] <- NA
  tao_imp <- impute_logreg(tao_imp, is_hot ~ sea_surface_temp)
  # Set humidity to NA in places where it was originally missing and re-impute it
  tao_imp$humidity[missing_humidity] <- NA
  tao_imp <- impute_lm(tao_imp, humidity ~ sea_surface_temp + air_temp)
}

Sesión 3

Imputación con bosques aleatorios

Un enfoque de aprendizaje automático para la imputación puede ser más preciso y más fácil de implementar en comparación con modelos estadísticos tradicionales. Primero, no requiere que especifiques relaciones entre variables. Además, los modelos de aprendizaje automático como los bosques aleatorios son capaces de descubrir relaciones altamente complejas y no lineales y explotarlas para predecir valores faltantes.

En este ejercicio, te familiarizarás con el paquete missForest, que construye un bosque aleatorio separado para predecir valores faltantes para cada variable, uno por uno. Llamarás a la función de imputación en los datos de películas biográficas, biopics, con los que has trabajado anteriormente en el curso y luego extraerás los datos completos, así como los errores de imputación estimados.

¡Plantemos algunos bosques aleatorios!

biopics <- read.csv("~/edi_imp/datacamp_JB/handing/biopics.csv")
# Load the missForest package
library(missForest)
## 
## Attaching package: 'missForest'
## The following object is masked from 'package:VIM':
## 
##     nrmse
# Impute biopics data using missForest
imp_res <- missForest(biopics)

# Extract imputed data and check for missing values
imp_data <- imp_res$ximp
print(sum(is.na(imp_data)))
## [1] 0
# Extract and print imputation errors
imp_err <- imp_res$OOBerror
print(imp_err)
##      NRMSE        PFC 
## 0.02041527 0.04033688

Variable-wise imputation errors

En el ejercicio anterior has extraído los errores de imputación estimados a partir de la salida de missForest. Esto te dio dos números:

el error cuadrático medio raíz normalizado (NRMSE) para todas las variables continuas; la proporción de entradas falsamente clasificadas (PFC) para todas las variables categóricas. Sin embargo, ¡podría darse el caso de que el modelo de imputación funcione muy bien para una variable continua y muy mal para otra! Para diagnosticar tales casos, basta con decirle a missForest que produzca estimaciones de error por variable. Esto se hace estableciendo el argumento variablewise en TRUE.

¡Los datos biopics y el paquete missForest ya han sido cargados para ti, así que echemos un vistazo más de cerca a los errores!

# Impute biopics data with missForest computing per-variable errors
imp_res <- missForest(biopics, variablewise = TRUE)

# Extract and print imputation errors
per_variable_errors <- imp_res$OOBerror
print(per_variable_errors)
##          PFC          MSE          MSE          MSE          PFC          PFC 
##    0.0000000    0.0000000 1309.5395901    0.0000000    0.0000000    0.1755319 
##          MSE          PFC 
##    0.0000000    0.0000000
# Rename errors' columns to include variable names
names(per_variable_errors) <- paste(names(biopics), 
                                    names(per_variable_errors),
                                    sep = "_")

# Print the renamed errors
print(per_variable_errors)
##   country_PFC      year_MSE  earnings_MSE   sub_num_MSE  sub_type_PFC 
##     0.0000000     0.0000000  1309.5395901     0.0000000     0.0000000 
##  sub_race_PFC non_white_MSE   sub_sex_PFC 
##     0.1755319     0.0000000     0.0000000

Observa cómo produjiste una serie de medidas de error en lugar de las dos por defecto que has visto antes. ¡Ahora puedes evaluar la calidad de imputación para cada variable por separado! Esto es útil cuando necesitas saber cómo se desempeña el modelo para una variable en particular que deseas modelar o analizar más a fondo.

Compromiso velocidad-precisión

En el último video, has visto que hay dos perillas que puedes ajustar para influir en el rendimiento de los bosques aleatorios:

Número de árboles de decisión en cada bosque. Número de variables utilizadas para la división dentro de los árboles de decisión. Aumentar cada uno de ellos puede mejorar la precisión del modelo de imputación, pero también requerirá más tiempo para ejecutarse. En este ejercicio, explorarás estas ideas por ti mismo ajustando missForest() a los datos de biopics dos veces con diferentes configuraciones. Mientras sigues las instrucciones, presta atención a los errores que imprimirás y al tiempo que tomará la ejecución del código.

# Set number of trees to 5 and number of variables used for splitting to 2
imp_res <- missForest(biopics, mtry = 2, ntree = 5)

# Print the resulting imputation errors
print(imp_res$OOBerror)
##      NRMSE        PFC 
## 0.02462457 0.08496094
# Set number of trees to 50 and number of variables used for splitting to 6
imp_res <- missForest(biopics, mtry = 6, ntree = 50)

# Print the resulting imputation errors
print(imp_res$OOBerror)
##      NRMSE        PFC 
## 0.02175938 0.04521277

Compara los errores y los tiempos de ejecución de los dos modelos de imputación. ¿Puedes ver una relación? Como dicen, “no hay nada gratuito”. Para obtener una imputación más precisa, tuviste que invertir más tiempo de computación. ¡Felicitaciones por terminar el capítulo! Nos vemos en el capítulo final, donde aprenderás a incorporar la incertidumbre de la imputación en tus análisis y predicciones.

Capítulo 4

Sesión 1

Envolver la imputación y el modelado en una función

Siempre que realice cualquier análisis o modelado en datos imputados, debe tener en cuenta la incertidumbre de la imputación. Ejecutar un modelo en un conjunto de datos imputados solo una vez ignora el hecho de que la imputación estima los valores faltantes con incertidumbre. Los errores estándar de dicho modelo tienden a ser demasiado pequeños. La solución a esto es la imputación múltiple y una forma de implementarla es mediante el bootstrap.

En los próximos ejercicios, trabajará con los datos familiares de biopics. El objetivo es utilizar la imputación múltiple mediante bootstrap y la regresión lineal para ver si, en función de los datos disponibles, las películas biográficas con mujeres ganan menos que las de hombres.

Comencemos escribiendo una función que cree una muestra de bootstrap, la impute y ajuste un modelo de regresión lineal.

calc_gender_coef <- function(data, indices) {
  # Get bootstrap sample
  data_boot <- data[indices, ]
  # Impute with kNN imputation
  data_imp <- kNN(data_boot, k = 5)
  # Fit linear regression
  linear_model <- lm(earnings ~ sub_sex + sub_type + year, data = data_imp)
  # Extract and return gender coefficient
  gender_coefficient <- coef(linear_model)[2]
  return(gender_coefficient)
}

La función calc_gender_coef() que acabas de programar toma los datos y los índices de bootstrap como entradas, y produce nuestra estadística de interés: el impacto del género en las ganancias de la regresión lineal. ¡Ahora puedes usar esta función en el algoritmo de bootstrapping!

Corriendo con bootstrap

Buen trabajo escribiendo calc_gender_coef() en el último ejercicio! Esta función crea una muestra de bootstrap, la imputa y produce el coeficiente de regresión lineal que describe el impacto de que el tema de la película sea femenino en las ganancias de la película.

En este ejercicio, usarás el paquete boot para obtener una distribución de bootstrap de estos coeficientes. La propagación de esta distribución capturará la incertidumbre de la imputación. También verás cómo la distribución de bootstrap difiere de una imputación y regresión de una sola vez. ¡Vamos a hacer un poco de bootstrapping!

# Load the boot library
library(boot)

# Run bootstrapping on biopics data
boot_results <- boot(biopics, statistic = calc_gender_coef, R = 50)

# Print and plot bootstrapping results
print(boot_results)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = biopics, statistic = calc_gender_coef, R = 50)
## 
## 
## Bootstrap Statistics :
##     original    bias    std. error
## t1* 2.060346 0.3033707    5.473554
plot(boot_results)

¡Excelente bootstrapping! Si hubieras ejecutado la imputación kNN y el análisis de regresión en los datos de biopics solo una vez, habrías obtenido un coeficiente de -1.45 para las películas sobre mujeres (llamado “original” en la salida de la consola), lo que sugiere que las películas sobre mujeres ganan menos. Sin embargo, al corregir la incertidumbre de la imputación, ¡has obtenido una distribución que cubre tanto valores negativos como positivos!

Bootstrapping intervalos de confianza

Después de haber generado la distribución del coeficiente del efecto femenino en el último ejercicio, ahora puedes usarla para estimar un intervalo de confianza. Esto te permitirá hacer la siguiente evaluación sobre tus datos: “Dada la incertidumbre de la imputación, estamos 95% seguros de que el efecto femenino en las ganancias se encuentra entre a y b”, donde a y b son los límites inferior y superior del intervalo.

En el último ejercicio, ejecutaste la técnica de bootstrapping con R = 50 réplicas. Sin embargo, en la mayoría de las aplicaciones esto no es suficiente. En este ejercicio, puedes utilizar los boot_results que se prepararon para ti utilizando 1000 réplicas. Primero, verás si la distribución de bootstrapping parece normal. Si es así, entonces podrás confiar en la distribución normal para calcular el intervalo de confianza.

# Run bootstrapping on biopics data
#boot_results <- boot(biopics, statistic = calc_gender_coef, R = 1000)
# Plot and print boot_results
plot(boot_results)

print(boot_results)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = biopics, statistic = calc_gender_coef, R = 50)
## 
## 
## Bootstrap Statistics :
##     original    bias    std. error
## t1* 2.060346 0.3033707    5.473554
# Calculate and print confidence interval
boot_ci <- boot.ci(boot_results, conf = 0.95, type = "norm")
print(boot_ci)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 50 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_results, conf = 0.95, type = "norm")
## 
## Intervals : 
## Level      Normal        
## 95%   (-8.971, 12.485 )  
## Calculations and Intervals on Original Scale

A pesar de que la tendencia general parece ser una relación negativa, las réplicas de bootstrap muestran que algunas películas con protagonistas femeninas en realidad ganan más. Al tener en cuenta la incertidumbre de la imputación, no se puede estar al 100% seguro acerca de la dirección de esta relación, aunque un análisis único sugiera lo contrario.

Sesión 2

El flujo de MICE: mice - with - pool

El flujo de MICE (imputación múltiple por ecuaciones encadenadas) nos permite estimar la incertidumbre de la imputación mediante la imputación de un conjunto de datos varias veces mediante la imputación basada en modelos, mientras se extrae de las distribuciones condicionales. De esta manera, cada conjunto de datos imputados es ligeramente diferente. Luego, se realiza un análisis en cada uno de ellos y se combinan los resultados, obteniendo las cantidades de interés junto con sus intervalos de confianza que reflejan la incertidumbre de la imputación.

En este ejercicio, practicarás el flujo típico de MICE: mice() - with() - pool(). Realizarás un análisis de regresión en los datos de biopics para ver qué tipo de ocupación de sujeto, sub_type, está asociada con mayores ingresos en películas. ¡Juguemos con MICE!

# Load mice package
library(mice)
## 
## Attaching package: 'mice'
## The following object is masked from 'package:stats':
## 
##     filter
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
# Impute biopics with mice using 5 imputations
biopics_multiimp <- mice(biopics, m = 5, seed = 3108)
## 
##  iter imp variable
##   1   1  earnings  sub_race
##   1   2  earnings  sub_race
##   1   3  earnings  sub_race
##   1   4  earnings  sub_race
##   1   5  earnings  sub_race
##   2   1  earnings  sub_race
##   2   2  earnings  sub_race
##   2   3  earnings  sub_race
##   2   4  earnings  sub_race
##   2   5  earnings  sub_race
##   3   1  earnings  sub_race
##   3   2  earnings  sub_race
##   3   3  earnings  sub_race
##   3   4  earnings  sub_race
##   3   5  earnings  sub_race
##   4   1  earnings  sub_race
##   4   2  earnings  sub_race
##   4   3  earnings  sub_race
##   4   4  earnings  sub_race
##   4   5  earnings  sub_race
##   5   1  earnings  sub_race
##   5   2  earnings  sub_race
##   5   3  earnings  sub_race
##   5   4  earnings  sub_race
##   5   5  earnings  sub_race
## Warning: Number of logged events: 50
# Fit linear regression to each imputed data set 
lm_multiimp <- with(biopics_multiimp, lm(earnings ~ year + sub_type))

# Pool and summarize regression results
lm_pooled <- pool(lm_multiimp)
summary(lm_pooled, conf.int = TRUE, conf.level = 0.95)
##                              term     estimate  std.error   statistic
## 1                     (Intercept) -287.8866750 490.062824 -0.58744851
## 2                            year    0.1612176   0.239277  0.67376974
## 3  sub_typeAcademic (Philosopher)  -32.8423364  39.593926 -0.82947915
## 4                sub_typeActivist  -17.4281787  16.052579 -1.08569339
## 5                   sub_typeActor  -30.7740287  19.378762 -1.58802862
## 6                 sub_typeActress  -27.3033456  21.249448 -1.28489670
## 7      sub_typeActress / activist   16.0962907  39.192849  0.41069458
## 8                  sub_typeArtist  -27.4779956  18.148136 -1.51409465
## 9                 sub_typeAthlete   -4.2336944  12.503822 -0.33859202
## 10     sub_typeAthlete / military   79.1943734  38.055447  2.08102595
## 11                 sub_typeAuthor  -23.6540827  18.471621 -1.28056347
## 12          sub_typeAuthor (poet)  -23.4506193  19.985477 -1.17338304
## 13               sub_typeComedian  -22.9330598  21.576033 -1.06289508
## 14               sub_typeCriminal   -2.6478096  16.754147 -0.15803906
## 15             sub_typeGovernment   -5.0221576  21.879188 -0.22954040
## 16             sub_typeHistorical   -9.3734433  19.183957 -0.48860845
## 17             sub_typeJournalist  -26.5071032  29.005345 -0.91386960
## 18                  sub_typeMedia  -11.5302020  26.461566 -0.43573393
## 19               sub_typeMedicine    4.7788761  15.006237  0.31845932
## 20               sub_typeMilitary   10.9417685  19.325322  0.56618816
## 21    sub_typeMilitary / activist   37.2248141  39.752579  0.93641257
## 22               sub_typeMusician  -19.7931797  17.320867 -1.14273611
## 23                  sub_typeOther  -15.5895605  16.044509 -0.97164460
## 24             sub_typePolitician  -13.6751859  39.752579 -0.34400752
## 25                 sub_typeSinger    0.6677979  17.844592  0.03742298
## 26                sub_typeTeacher   51.1575084  39.268925  1.30274786
## 27           sub_typeWorld leader    5.9976436  16.882980  0.35524793
##            df    p.value        2.5 %       97.5 %
## 1    3.983993 0.58858265 -1650.677928 1074.9045785
## 2    4.030421 0.53712441    -0.501149    0.8235843
## 3  139.785256 0.40824758  -111.122704   45.4380311
## 4   10.576004 0.30174143   -52.933253   18.0768961
## 5   10.847500 0.14097812   -73.499669   11.9516115
## 6    7.876739 0.23532207   -76.438456   21.8317642
## 7  169.869958 0.68181408   -61.271473   93.4640548
## 8    8.280382 0.16719720   -69.082290   14.1262992
## 9   15.246283 0.73953513   -30.847513   22.3801244
## 10 336.810947 0.03818668     4.338081  154.0506662
## 11   7.047583 0.24087818   -67.272814   19.9646489
## 12  10.612273 0.26630514   -67.635233   20.7339945
## 13  16.686775 0.30297041   -68.519689   22.6535693
## 14   7.275709 0.87872330   -41.962758   36.6671385
## 15  81.573759 0.81902349   -48.550237   38.5059216
## 16   6.172715 0.64199292   -55.998343   37.2514567
## 17 113.284343 0.36272632   -83.970362   30.9561559
## 18   6.128774 0.67795894   -75.950965   52.8905609
## 19 247.938063 0.75040464   -24.777080   34.3348320
## 20   6.645835 0.58986247   -35.253693   57.1372305
## 21 130.043524 0.35079655   -41.420661  115.8702894
## 22   6.996299 0.29074124   -60.754913   21.1685537
## 23   7.168854 0.36286255   -53.348394   22.1692725
## 24 130.043524 0.73139625   -92.320661   64.9702894
## 25  10.377953 0.97085786   -38.897026   40.2326216
## 26 163.415623 0.19449369   -26.382404  128.6974204
## 27  14.000056 0.72769899   -30.212733   42.2080205

En este caso, has seguido el flujo “mice-with-pool” para imputar, modelar y agrupar los resultados. Ahora, echa un vistazo a la salida en la consola: algunos sub_types tienen un impacto positivo en las ganancias. Sin embargo, al tener en cuenta la incertidumbre de la imputación con una confianza del 95%, nunca estamos seguros de estos efectos, ¡ya que los límites inferiores son negativos! Con una excepción: para sub_typeAthlete / military, tanto los límites inferiores como los superiores son positivos. Lo que podemos decir con seguridad es que ¡las películas sobre atletas militares son populares!

Selección de modelos por defecto

MICE crea un modelo de imputación separado para cada variable en los datos. El tipo de modelo depende del tipo de variable en cuestión. Una forma popular de especificar los tipos de modelos que queremos usar es establecer un modelo predeterminado para cada uno de los cuatro tipos de variables.

Puede hacer esto pasando el argumento defaultMethod a mice(), que debe ser un vector de longitud 4 que contenga los métodos de imputación predeterminados para:

Variables continuas, Variables binarias, Variables categóricas (factores no ordenados), Variables factoriales (factores ordenados).

En este ejercicio, aprovechará la documentación de mice para ver la lista de métodos disponibles y seleccionar los deseados para que el algoritmo los use. ¡Vamos a hacer selección de modelos!

# Impute biopics using the methods specified in the instruction
biopics_multiimp <- mice(biopics, m = 20, 
                         defaultMethod = c("cart", "lda", "pmm", "polr"))
## 
##  iter imp variable
##   1   1  earnings  sub_race
##   1   2  earnings  sub_race
##   1   3  earnings  sub_race
##   1   4  earnings  sub_race
##   1   5  earnings  sub_race
##   1   6  earnings  sub_race
##   1   7  earnings  sub_race
##   1   8  earnings  sub_race
##   1   9  earnings  sub_race
##   1   10  earnings  sub_race
##   1   11  earnings  sub_race
##   1   12  earnings  sub_race
##   1   13  earnings  sub_race
##   1   14  earnings  sub_race
##   1   15  earnings  sub_race
##   1   16  earnings  sub_race
##   1   17  earnings  sub_race
##   1   18  earnings  sub_race
##   1   19  earnings  sub_race
##   1   20  earnings  sub_race
##   2   1  earnings  sub_race
##   2   2  earnings  sub_race
##   2   3  earnings  sub_race
##   2   4  earnings  sub_race
##   2   5  earnings  sub_race
##   2   6  earnings  sub_race
##   2   7  earnings  sub_race
##   2   8  earnings  sub_race
##   2   9  earnings  sub_race
##   2   10  earnings  sub_race
##   2   11  earnings  sub_race
##   2   12  earnings  sub_race
##   2   13  earnings  sub_race
##   2   14  earnings  sub_race
##   2   15  earnings  sub_race
##   2   16  earnings  sub_race
##   2   17  earnings  sub_race
##   2   18  earnings  sub_race
##   2   19  earnings  sub_race
##   2   20  earnings  sub_race
##   3   1  earnings  sub_race
##   3   2  earnings  sub_race
##   3   3  earnings  sub_race
##   3   4  earnings  sub_race
##   3   5  earnings  sub_race
##   3   6  earnings  sub_race
##   3   7  earnings  sub_race
##   3   8  earnings  sub_race
##   3   9  earnings  sub_race
##   3   10  earnings  sub_race
##   3   11  earnings  sub_race
##   3   12  earnings  sub_race
##   3   13  earnings  sub_race
##   3   14  earnings  sub_race
##   3   15  earnings  sub_race
##   3   16  earnings  sub_race
##   3   17  earnings  sub_race
##   3   18  earnings  sub_race
##   3   19  earnings  sub_race
##   3   20  earnings  sub_race
##   4   1  earnings  sub_race
##   4   2  earnings  sub_race
##   4   3  earnings  sub_race
##   4   4  earnings  sub_race
##   4   5  earnings  sub_race
##   4   6  earnings  sub_race
##   4   7  earnings  sub_race
##   4   8  earnings  sub_race
##   4   9  earnings  sub_race
##   4   10  earnings  sub_race
##   4   11  earnings  sub_race
##   4   12  earnings  sub_race
##   4   13  earnings  sub_race
##   4   14  earnings  sub_race
##   4   15  earnings  sub_race
##   4   16  earnings  sub_race
##   4   17  earnings  sub_race
##   4   18  earnings  sub_race
##   4   19  earnings  sub_race
##   4   20  earnings  sub_race
##   5   1  earnings  sub_race
##   5   2  earnings  sub_race
##   5   3  earnings  sub_race
##   5   4  earnings  sub_race
##   5   5  earnings  sub_race
##   5   6  earnings  sub_race
##   5   7  earnings  sub_race
##   5   8  earnings  sub_race
##   5   9  earnings  sub_race
##   5   10  earnings  sub_race
##   5   11  earnings  sub_race
##   5   12  earnings  sub_race
##   5   13  earnings  sub_race
##   5   14  earnings  sub_race
##   5   15  earnings  sub_race
##   5   16  earnings  sub_race
##   5   17  earnings  sub_race
##   5   18  earnings  sub_race
##   5   19  earnings  sub_race
##   5   20  earnings  sub_race
## Warning: Number of logged events: 200
# Print biopics_multiimp
print(biopics_multiimp)
## Class: mids
## Number of multiple imputations:  20 
## Imputation methods:
##   country      year  earnings   sub_num  sub_type  sub_race non_white   sub_sex 
##        ""        ""    "cart"        ""        ""     "pmm"        ""        "" 
## PredictorMatrix:
##          country year earnings sub_num sub_type sub_race non_white sub_sex
## country        0    1        1       1        1        1         1       1
## year           1    0        1       1        1        1         1       1
## earnings       1    1        0       1        1        1         1       1
## sub_num        1    1        1       0        1        1         1       1
## sub_type       1    1        1       1        0        1         1       1
## sub_race       1    1        1       1        1        0         1       1
## Number of logged events:  200 
##   it im      dep meth                            out
## 1  1  1 earnings cart sub_typeAcademic (Philosopher)
## 2  1  1 sub_race  pmm    sub_typeMilitary / activist
## 3  1  2 earnings cart sub_typeAcademic (Philosopher)
## 4  1  2 sub_race  pmm    sub_typeMilitary / activist
## 5  1  3 earnings cart sub_typeAcademic (Philosopher)
## 6  1  3 sub_race  pmm    sub_typeMilitary / activist

La capacidad de especificar modelos de imputación puede resultar útil cuando se observa que algunos métodos específicos no funcionan bien. Otro factor que influye en cómo funcionan los métodos de imputación es el conjunto de predictores que utilizan. En el siguiente ejercicio, veremos cómo establecer estos predictores.

Usando una matriz predictora

Se trata de tomar decisiones importantes cuando se utiliza la imputación basada en modelos, como por ejemplo, qué variables deben incluirse como predictores y en qué modelos. En mice(), esto se rige por la matriz de predictores y, por defecto, todas las variables se utilizan para imputar todas las demás.

En caso de tener muchas variables en los datos o poco tiempo para realizar una selección adecuada del modelo, puede utilizar la funcionalidad de mice para crear una matriz de predictores basada en las correlaciones entre las variables. Esta matriz se puede pasar entonces a mice(). En este ejercicio, practicará exactamente esto: primero construirá una matriz de predictores de modo que cada variable se impute utilizando las variables más correlacionadas con ella; luego, alimentará su matriz de predictores a la función de imputación. ¡Vamos a intentar esta simple selección de modelos!

# Create predictor matrix with minimum correlation of 0.1
pred_mat <- quickpred(biopics, mincor = 0.1)

# Impute biopics with mice
biopics_multiimp <- mice(biopics, 
                         m = 10, 
                         predictorMatrix = pred_mat,
                         seed = 3108)
## 
##  iter imp variable
##   1   1  earnings  sub_race
##   1   2  earnings  sub_race
##   1   3  earnings  sub_race
##   1   4  earnings  sub_race
##   1   5  earnings  sub_race
##   1   6  earnings  sub_race
##   1   7  earnings  sub_race
##   1   8  earnings  sub_race
##   1   9  earnings  sub_race
##   1   10  earnings  sub_race
##   2   1  earnings  sub_race
##   2   2  earnings  sub_race
##   2   3  earnings  sub_race
##   2   4  earnings  sub_race
##   2   5  earnings  sub_race
##   2   6  earnings  sub_race
##   2   7  earnings  sub_race
##   2   8  earnings  sub_race
##   2   9  earnings  sub_race
##   2   10  earnings  sub_race
##   3   1  earnings  sub_race
##   3   2  earnings  sub_race
##   3   3  earnings  sub_race
##   3   4  earnings  sub_race
##   3   5  earnings  sub_race
##   3   6  earnings  sub_race
##   3   7  earnings  sub_race
##   3   8  earnings  sub_race
##   3   9  earnings  sub_race
##   3   10  earnings  sub_race
##   4   1  earnings  sub_race
##   4   2  earnings  sub_race
##   4   3  earnings  sub_race
##   4   4  earnings  sub_race
##   4   5  earnings  sub_race
##   4   6  earnings  sub_race
##   4   7  earnings  sub_race
##   4   8  earnings  sub_race
##   4   9  earnings  sub_race
##   4   10  earnings  sub_race
##   5   1  earnings  sub_race
##   5   2  earnings  sub_race
##   5   3  earnings  sub_race
##   5   4  earnings  sub_race
##   5   5  earnings  sub_race
##   5   6  earnings  sub_race
##   5   7  earnings  sub_race
##   5   8  earnings  sub_race
##   5   9  earnings  sub_race
##   5   10  earnings  sub_race
## Warning: Number of logged events: 50
# Print biopics_multiimp
print(biopics_multiimp)
## Class: mids
## Number of multiple imputations:  10 
## Imputation methods:
##   country      year  earnings   sub_num  sub_type  sub_race non_white   sub_sex 
##        ""        ""     "pmm"        ""        "" "polyreg"        ""        "" 
## PredictorMatrix:
##          country year earnings sub_num sub_type sub_race non_white sub_sex
## country        0    0        0       0        0        0         0       0
## year           0    0        0       0        0        0         0       0
## earnings       1    1        0       0        0        1         1       0
## sub_num        0    0        0       0        0        0         0       0
## sub_type       0    0        0       0        0        0         0       0
## sub_race       0    1        1       1        1        0         1       0
## Number of logged events:  50 
##   it im      dep    meth                         out
## 1  1  1 sub_race polyreg sub_typeMilitary / activist
## 2  1  2 sub_race polyreg sub_typeMilitary / activist
## 3  1  3 sub_race polyreg sub_typeMilitary / activist
## 4  1  4 sub_race polyreg sub_typeMilitary / activist
## 5  1  5 sub_race polyreg sub_typeMilitary / activist
## 6  1  6 sub_race polyreg sub_typeMilitary / activist

Sesión 3

Analizando los patrones de datos faltantes

El primer paso para trabajar con datos incompletos es obtener información sobre los patrones de ausencia de datos, y una buena manera de hacerlo es mediante visualizaciones. Comenzarás tu análisis de los datos de África empleando el paquete VIM para crear dos visualizaciones: el gráfico de agregación y el gráfico de columna vertebral. Te dirán cuántos datos faltan, en qué variables y configuraciones, y si podemos decir algo sobre el mecanismo de ausencia de datos. ¡Comencemos con algunas gráficas!

africa <- read.csv("~/edi_imp/datacamp_JB/handing/africa.csv", sep = ";")
# Load VIM
library(VIM)



# Draw a combined aggregation plot of africa
africa %>%
  aggr(combined = TRUE, numbers = TRUE)

# Draw a spine plot of country vs trade
africa %>% 
  select(country, trade) %>%
  spineMiss()

Pregunta

Basándose en el gráfico de columna vertebral que acaba de crear, ¿cuál de las siguientes afirmaciones es VERDADERA?

Posibles respuestas

  1. Hay más datos faltantes en el comercio para Camerún que para Burundi.

  2. Hay más datos faltantes en gdp_pc para Burundi que para Camerún.

  3. No hay ningún país con más del 20% de valores faltantes en comercio.

Correcto, ¡no hay tantos valores faltantes! Además, observe en el gráfico de columna vertebral que los datos de África parecen ser MAR - al menos con respecto al PIB y al país, lo que significa que se pueden imputar.

###Imputando e inspeccionando resultados

¡Buen trabajo al visualizar los datos faltantes en el ejercicio anterior! Ha descubierto que hay algunas entradas faltantes en el PIB, gdp_pc, y en el comercio como porcentaje del PIB, comercio. Además, sospecha que los datos son MAR, por lo que se pueden imputar. En este ejercicio, hará uso de la imputación múltiple del paquete mice para imputar los datos de africa. Luego, dibujará un gráfico de tiras degdp_pc vs trade para ver si los datos imputados no rompen la relación entre estas variables. ¡Deje que mice haga el trabajo!

# Load mice
library(mice)

# Impute africa with mice
africa_multiimp <- mice(africa, m = 5, defaultMethod = "cart", seed = 3108)
## 
##  iter imp variable
##   1   1  gdp_pc  trade
##   1   2  gdp_pc  trade
##   1   3  gdp_pc  trade
##   1   4  gdp_pc  trade
##   1   5  gdp_pc  trade
##   2   1  gdp_pc  trade
##   2   2  gdp_pc  trade
##   2   3  gdp_pc  trade
##   2   4  gdp_pc  trade
##   2   5  gdp_pc  trade
##   3   1  gdp_pc  trade
##   3   2  gdp_pc  trade
##   3   3  gdp_pc  trade
##   3   4  gdp_pc  trade
##   3   5  gdp_pc  trade
##   4   1  gdp_pc  trade
##   4   2  gdp_pc  trade
##   4   3  gdp_pc  trade
##   4   4  gdp_pc  trade
##   4   5  gdp_pc  trade
##   5   1  gdp_pc  trade
##   5   2  gdp_pc  trade
##   5   3  gdp_pc  trade
##   5   4  gdp_pc  trade
##   5   5  gdp_pc  trade
# Draw a stripplot of gdp_pc versus trade
stripplot(africa_multiimp, gdp_pc ~ trade | .imp, pch = 20, cex = 1)

Se observa que la imputación funciona bien: hay pequeños grupos en los gráficos de dispersión, que probablemente corresponden a diferentes países. Cada punto de datos imputado encaja en uno de los grupos, en lugar de ser un valor atípico en algún lugar entre los grupos. Después de haber realizado la imputación, ¡ahora puedes proceder con el modelado!

Inferencia con datos imputados

En el último ejercicio, has utilizado mice para imputar los datos de africa. En este, implementarás los otros dos pasos del flujo de “mice - with - pool” que has aprendido anteriormente en el curso. El modelo de interés es una regresión lineal que explica el PIB, gdp_pc, con otras variables. Te interesa particularmente el coeficiente de libertades civiles, civlib. ¿Está más libertad asociada con un mayor crecimiento económico una vez que incorporamos la incertidumbre de la imputación? ¡Descubrámoslo!

# Fit linear regression to each imputed data set
lm_multiimp <- with(africa_multiimp, lm(gdp_pc ~ country + year + trade + infl + civlib))

# Pool regression results
lm_pooled <- pool(lm_multiimp)

# Summarize pooled results
summary(lm_pooled, conf.int = TRUE, conf.level = 0.9)
##               term      estimate    std.error  statistic       df      p.value
## 1      (Intercept) -30005.113539 5592.2088438 -5.3655209 107.8640 4.654452e-07
## 2   countryBurundi     70.655824   61.9717253  1.1401300 107.1336 2.567747e-01
## 3  countryCameroon    653.969299   58.3491217 11.2078688 106.1787 0.000000e+00
## 4     countryCongo   1337.668745  112.9974527 11.8380434 107.8509 0.000000e+00
## 5   countrySenegal    522.377263   78.6118111  6.6450226 107.4387 1.293220e-09
## 6    countryZambia    410.935079   83.5292689  4.9196537 107.8369 3.128698e-06
## 7             year     15.299939    2.8128543  5.4392932 107.8559 3.368002e-07
## 8            trade      5.059562    1.5814802  3.1992573 107.4800 1.811115e-03
## 9             infl     -4.347048    0.9837546 -4.4188336 108.0143 2.369012e-05
## 10          civlib    -40.661539  132.1221617 -0.3077571 106.5790 7.588678e-01
##              5 %          95 %
## 1  -39283.165341 -20727.061736
## 2     -32.167739    173.479388
## 3     557.148769    750.789829
## 4    1150.194108   1525.143382
## 5     391.947680    652.806846
## 6     272.351096    549.519062
## 7      10.633121     19.966758
## 8       2.435642      7.683482
## 9      -5.979179     -2.714917
## 10   -259.888744    178.565667

Pregunta

Basándose en el resumen de los resultados de la regresión agrupada que acaba de imprimir en la consola, ¿cuál de las siguientes afirmaciones sobre las libertades civiles en África es falsa?

Respuestas posibles

  1. En promedio, más libertad está asociada con un PIB más alto.

  2. Hay un 5% de probabilidad de que el coeficiente de civlib en la regresión de gdp_pc sea mayor que 342.154078.

  3. Basado en el intervalo de confianza del 90%, estamos seguros de que el impacto de civlib en gdp_pc es positivo.

Correcto, ¡esta es falsa! Dado que los límites inferior y superior tienen signos diferentes, no podemos estar seguros de la dirección del efecto. Felicitaciones, has recorrido un largo camino y aprendido mucho. ¡Bien hecho! Resumamos todo en el último video del curso.

Comentarios finales

  1. Observaciones finales ¡Felicidades por llegar hasta el final del curso! Recapitulemos lo que has aprendido.

  2. Lo que sabes En el Capítulo 1, viste cómo modelar datos incompletos puede ser problemático y que requiere un tratamiento especial. También aprendiste acerca de los tres mecanismos de datos faltantes y cómo obtener información sobre los patrones de datos faltantes mediante visualizaciones del paquete VIM y pruebas estadísticas. En el Capítulo 2, cubrimos la imputación basada en donantes. Primero, viste por qué la imputación media es típicamente una mala elección. Luego, aprendiste a utilizar la imputación de hot-deck y kNN del paquete VIM, junto con algunos trucos que hacen que funcionen aún mejor.

  3. Lo que sabes En el Capítulo 3, aprendiste el enfoque de imputación basado en modelos de hacer un bucle sobre variables e imputarlas hasta la convergencia. También viste cómo aumentar la variabilidad de los datos imputados dibujando de distribuciones condicionales. Finalmente, aprendiste acerca de la imputación basada en árboles con el paquete missForest. En el Capítulo 4, viste dos métodos para incorporar la incertidumbre de la imputación en el modelado: el bootstrapping utilizando el paquete boot y la imputación múltiple mediante ecuaciones encadenadas utilizando el paquete mice.

  4. ¿Qué método de imputación elegir? Aprendiste sobre muchos métodos diferentes de imputación. ¿Cuál elegir y cuándo? Aquí hay algunas pautas generales. Si tienes muchos datos o si tu imputación debe ejecutarse en tiempo real en producción, es mejor usar hot-deck rápido. Si sospechas relaciones específicas entre las variables basadas en el conocimiento del dominio, puedes usar este conocimiento en la imputación basada en modelos. De lo contrario, si la imputación no necesita ser muy rápida y las relaciones entre las variables no son obvias, un enfoque de aprendizaje automático como kNN o imputación basada en árboles es tu mejor opción.

  5. ¿Cómo estimar la incertidumbre de la imputación? También aprendiste sobre dos métodos para estimar la incertidumbre de la imputación: bootstrapping y mice. ¿Cuál elegir? Una vez más, permíteme ofrecer algunas pautas generales. Si tu aplicación debe ser relativamente rápida o si tienes ideas sobre qué modelos usar y cómo especificarlos, entonces mice es la mejor opción. De lo contrario, si deseas utilizar un método no paramétrico como kNN o hot-deck, o simplemente no quieres preocuparte por los supuestos de modelos específicos, el bootstrap podría ser una mejor opción.

  6. Siguientes pasos Una última cosa antes de que te vayas. Si deseas profundizar aún más en el conocimiento de la imputación, te sugiero que busques en Google “miceVignettes”. Los autores del paquete “mice” proporcionan una serie de seis viñetas con código R y ejemplos que abordan temas como la imputación pasiva, el post-procesamiento de datos imputados, la imputación de datos multinivel o el análisis de sensibilidad. Si absorbes fácilmente el conocimiento de los libros, entonces “Flexible Imputation of Missing Data” de Stef van Buuren es una lectura obligatoria. También utiliza el paquete “mice”. Finalmente, hay otros excelentes paquetes de R que vale la pena explorar y que no tuvimos tiempo de cubrir en este curso, como Amelia o mi, que permiten la imputación de series de tiempo y datos de panel.

  7. ¡Felicidades y buena suerte! Una vez más, felicidades por terminar el curso y gracias por acompañarme. Espero que el conocimiento y las habilidades que has adquirido te faciliten el trabajo con datos incompletos y te hagan más productivo. ¡Buena suerte!